In [1]:

    

%pylab inline
from sympy.interactive import init_printing
init_printing()
from sympy import pi, var, S, Piecewise, piecewise_fold
var("r R")
Vh = Piecewise((-S(2)/3 * pi * (3*R**2 - r**2), r <= R), (-S(4)/3 * pi * R**3 / r, True))
def laplace(f):
    return (r*f).diff(r, 2)/r
print "Vh ="
Vh


    

    




Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.
Vh =






    
Out[1]:





$$\begin{cases} - \frac{2}{3} \pi \left(3 R^{2} - r^{2}\right) & \text{for}\: r \leq R \\- \frac{4}{3} \frac{\pi R^{3}}{r} & \text{otherwise} \end{cases}$$

In [2]:

    

piecewise_fold(-laplace(Vh)/(4*pi))


    

    
Out[2]:





$$\begin{cases} -1 & \text{for}\: r \leq R \\0 & \text{otherwise} \end{cases}$$

In [3]:

    

from numpy import (empty, pi, meshgrid, linspace, sum, sin, exp, shape, sqrt,
        conjugate)
from numpy.fft import fftn, fftfreq, ifftn

N = 100
print "N =", N
L = 2.4*4
R = 1.
x1d = linspace(-L/2, L/2, N+1)[:-1]
x, y, z = meshgrid(x1d, x1d, x1d, indexing="ij")

r = sqrt(x**2+y**2+z**2)
nr = empty(shape(x), dtype="double")
nr[:] = 0
nr[r <= R] = -1

Vanalytic = empty(shape(x), dtype="double")
Vanalytic[r <= R] = -2./3 * pi * (3*R**2 - r[r <= R]**2)
Vanalytic[r > R] = -4./3 * pi * R**3 / r[r > R]

ng = fftn(nr) / N**3

G1d = N * fftfreq(N) * 2*pi/L
kx, ky, kz = meshgrid(G1d, G1d, G1d)
G2 = kx**2+ky**2+kz**2
G2[0, 0, 0] = 1  # omit the G=0 term

tmp = 2*pi*abs(ng)**2 / G2
tmp[0, 0, 0] = 0  # omit the G=0 term
E = sum(tmp) * L**3
print "Hartree Energy (calculated): %.15f" % E


Vg = 4*pi*ng / G2
Vg[0, 0, 0] = 0  # omit the G=0 term
V = ifftn(Vg).real * N**3
V += Vanalytic[N/2, N/2, N/2] - V[N/2, N/2, N/2]
l2_norm = sum((Vanalytic - V)**2) 
print "l2_norm = ", l2_norm
plot(x[:, N/2, N/2], Vanalytic[:, N/2, N/2], label="analytic")
plot(x[:, N/2, N/2], V[:, N/2, N/2], label="FFT")
legend(loc="best");


    

    




N = 100
Hartree Energy (calculated): 7.945878466121926
l2_norm =  48304.2303309






    

In [3]: